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Spurious velocities are unphysical currents that appear close to curved interfaces in diffuse inter- 
face methods. We analyse the causes of these spurious velocities in the free energy lattice Boltzmann 
algorithm. By making a suitable choice of the equilibrium distribution, and by finding the best way 
to numerically calculate derivatives, we show that these velocities may be decreased by an order of 
magnitude compared to previous models. Furthermore, we propose a momentum conserving forcing 
method that reduces spurious velocities by another factor of ~ 5. In three dimensions we find that 
19 velocity vectors is the minimum number necessary. 



I. INTRODUCTION 

A commonly used approach for the simulation of multi- 
phase fluid dynamics is the free energy lattice Boltzmann 
method introduced by Swift et al. This constitutes a 
so-called mesoscale method because it numerically solves 
the continuum equations of fluid dynamics by exploiting 
the underlying microscopic structure of these equations, 
without resorting to a description of the fluid in terms 
of molecular dynamics. One obstacle to simulating some 
systems is that discretisation errors lead to unphysical 
flows near interfaces. These so-called spurious velocities 
are present in multi-phase lattice Boltzmann methods 
and in the other diffuse interface methods. 

An illustration of these spurious velocities is given in 
Fig. [ifa), which shows the flow profile around a liquid 
drop coexisting with a surrounding gas phase. The sim- 
ulation is left until the long time steady state behaviour 
is reached. From a physical point of view all velocities 
should go to zero. What is observed, however, is that 
spurious flows persist indeflnitely. 

A number of papers have dealt with this problem. 
Wagner 0] analysed the case of binary fluids, and iden- 
tified that one way to eradicate spurious velocities was 
to remove non-ideal terms from the pressure tensor and 
introduce these as a body force of the form g = — t/fV/x^. 
However, because this is no longer written in terms of 
the divergence of a pressure tensor (note that in general 
— 0V/i^ can always be rewritten as —dpPap) then mo- 
mentum is no longer conserved. Furthermore, Wagner 
pointed out that this method is numerically unstable un- 
less some additional viscosity is artificially added to the 
system. 

Lee and Fisher [3] use another forcing method for a 
different implementation of the lattice Boltzmann algo- 
rithm. Again, they eliminate spurious velocities at the 
expense of sacrificing momentum conservation. An addi- 
tional difficulty with using forcing methods (including the 
one we present later) is that in order to update each lat- 
tice site the algorithm requires information from 2 lattice 
sites away, rather than just 1 for the standard method. 
This makes boundary conditions more complicated and 
slows down parallel computations, since more informa- 
tion needs to be passed between processors. 



Seta and Okui [4|] used a lattice Boltzmann scheme pro- 
posed by Inamuro et al. and considered calculating 
the derivatives in the pressure tensor using a more accu- 
rate fourth-order scheme (as opposed to the usual second 
order accurate method). As will be shown later, however, 
they do not choose an optimum equilibrium distribution 
and hence their improvement in the spurious velocities is 
limited. 

In this paper, we analyse the free energy lattice Boltz- 
mann scheme for a liquid-gas system proposed by Swift 
et al. [l| and show that by making a careful choice of 
the equilibrium distribution (and also finding the best 
way to calculate derivatives) the magnitude of spurious 
velocities can be significantly reduced. Furthermore, we 
present a second numerical scheme which moves gradient 
terms in the equilibrium distribution into a body force. 
This leads to a further reduction in spurious velocities 
whilst preserving momentum conservation. 

The results of this analysis are equally applicable to 
other multiphase systems {e.g. binary fiuids) when the 
free energy lattice Boltzmann method is used to solve 
their equations of motion. 

II. MODEL 

The pressure tensor for a liquid-gas system using a 
Landau free energy is given by 

PaP = [PO - KpV^p - -^IV/Op) Sap + ndapdpp, (1) 

where p is the fluid density and k is a parameter related 
to the surface tension. We choose the bulk pressure pq 
to be that of a van der Waals fluid, 

PT 2 fo\ 

P' = T^p~''P- 

This leads to liquid-gas phase separation below a critical 
temperature. 

The analysis which follows is performed for a D2Q9 
lattice Boltzmann scheme (in section IVIII the results are 
summarised for the D3Q19 model) which uses a square 
lattice of side Ax, time-step At, and has 9 velocity vec- 
tors, Bi, where Gq = (0,0), ei,2 = (±c, 0), 63,4 = (0,±c), 
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e5,6 = (±c, ±c), and ey^g = (=Fc, ±c). The parameter 
c = is a lattice velocity. 

A particle distribution function fi{r,t) gives the mass 
density of particles travelling from lattice site r, at time 
i, in a direction e^. The physical variables are related to 
this distribution function by 



pu. 



(3) 



where p is the mass density and u is the velocity of the 
fluid. 

The time evolution equation for the particle distribu- 
tion function, using the standard BGK approximation, is 
given by 

Mr + e,Ai, t + At) = /,(r, t) - ^ [/, - + F,, (4) 

where r is a relaxation parameter related to the viscos- 
ity, and f^"^ is an equilibrium distribution. It has been 
shown previously that this reduces to the Navier-Stokes 
equation provided the moments of f^'' and Fi are chosen 
suitably 6] (see appendix[A| . The final Fi term is respon- 
sible for introducing a body force. This is not present in 
the standard formulation of the free energy lattice Boltz- 
mann algorithm and so for now we set it to zero. In 
section I VII however, we discuss how this term can be 
usefully implemented to help reduce spurious velocities 
further. 

The equilibrium distribution can be written as 



(3 r £■ 

[pUaUp + A [Uad/3P + UpdaP + Sa^U^fd-^p]) 

+ ji ( w'fpo - w*pV^p + wf'K.d^pd^p 

+ ndypdyp^ w"!^ ndxpdypj, (5) 



for i = 1, .., 8, where u'i_4 = i, w^.s, = tk' ^'^'^ summation 
over repeated indices is assumed. The i = stationary 
value is chosen to conserve mass: 



r,\v) = p-Y,u{v). 



(6) 



The top two lines on the right hand side of Eq. ([5|) cor- 
respond to a standard expansion of the Maxwell Boltz- 
mann distribution in discretised space 0], and a cor- 
rection term involving A (see Eq. (jA4p ) which ensures 
Galilean invariance 6]. These terms are not important 
from the point of view of spurious velocities because they 
each contain the fluid velocity to some power, which 
is expected to be zero in equilibrium. 

The last two lines in Eq. ([5]) give the pressure tensor 
contribution to the equilibrium distribution. This has 
been written in its most general form involving the free 
parameter weights , w\, wf^, wf^, and . Through 



the course of this paper optimum values for these param- 
eters will be obtained. 

The derivatives in the equilibrium distribution ^ are 
explicitly calculated within the algorithm using finite dif- 
ference schemes. For instance, one simple choice for cal- 
culating the X derivative of p is given by 



(7) 



The bar above the partial derivative denotes that this is 
a discrete operator. By Taylor expanding the right hand 
side we find that 



(8) 



The discrete operator is correct up to second order but 
there are higher order terms which are responsible for 
generating the spurious flows. 

A useful representation of finite difference operators is 
to denote them by stencils. For instance Eq. ([7]) can be 
rewritten 



dxP 



1 



2Ax 












1 




(9) 



The central entry in the matrix represents the point at 
which the derivative is being made and the surrounding 
8 entries correspond to the neighbouring lattice points 
surrounding this. This, however, is not the only choice 
for calculating the x derivative. The most general stencil 
using only 9 lattice nodes can be written 



dx = 



1 
Ax 

a. 4 



-BOB 
-A A 
-BOB 



^Ax^dl + 2BAx^d^d^ + 



(10) 



where i? is a free parameter which can be used to deter- 
mine the third order term and A is defined by 2 A + AB — 
1. 

Similarly, the Laplacian operator can be represented 

by 



1 



Ax" 



D C D 

C -A{C + D) C 
D C D 



p 



= V^ + ^{dt + d^)+DAx^dld.^ 



.(11) 



where C + 2D=l. 

In equilibrium the Navier-Stokes equation reduces to 







{12) 



In terms of the lattice Boltzmann algorithm, the partial 
derivative operator acting on the pressure tensor in Eq. 

is implemented as a result of the choice of equilib- 
rium distribution and the streaming and colliding oper- 
ations. When T — 1 (in section |V] we discuss the more 
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general case) the lattice Boltzmann equation ^ reduces 
to 



/,(r,i + Ai) = /f'(r-e,At,i) 



(13) 



We consider the idealised case when at some time t 
the system is at rest, i.e. u(r,t) — 0, and the density 
distribution is chosen such that the continuous opera- 
tor equation (fT2|) is solved exactly. We ask the ques- 
tion what happens when the continuous operators are 
replaced by their discrete counterparts. In this case 
will no longer be exactly satisfied and instead there will 
be some spurious force G on the left hand side. 

Using Eq. (O, this force can be expressed in terms of 
stencils of the various terms in the equilibrium distribu- 
tion: 



^ t + At)~ pu^ir, <)) 



w 



These constraints are also 
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5-6 ^ 4' and Wj_% = -J. 
necessary to obtain the correct moments of the equilib- 
rium distribution in Eq. (jA2|. 

Given all these conditions the spurious force can be 
rewritten as 
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in, for example, the x direction. Here, we define Map = 
ndapdpp. In writing this we have made use of the symme- 
try properties of the system to immediately reduce the 
number of free parameters in the model. For instance, 
the bulk pressure po does not have a preferred direction 
[i.e. it acts the same in the x and y directions) and hence 



we expect that = 

V 



t(;J_4, and 



P 



w: 



W4, which we denote by 



VJ 



P 



Other terms do 



XX 

A 



have a preferred direction. For example, M^x is less re- 
stricted and has the constraints wf^ = w^^, w'^^ — w, 
and wf^ = Wq^ = w^^ = Wg'^. Because the equilibrium 
should be invariant under simultaneous interchange of x 
and y and switching the velocities 61^2,7 ©3,4, 8) then we 

mXX — nn^y aud VJc:"^ — 



"1-2' 



expect that Wi^2 = '^l-i^ ^3-4 = "^i 

To first order, Gx should agree with Eq. 
in the x direction is given by 



5-8- 



dni), which 



Gx 



-dx {po-KpV'^p)-^dx (Mxx- Myy) - dyMxy.{15) 



By comparing Eq. (|15p with Eq. further restric- 

tions are possible. For instance, by using Eq. pO)) . 
the Po stencil becomes —dxPo to second order provided 
that 2u;^4 -\- Aw^ ^ = 1. Similarly, 2w{ 4 + 4wl g = 1, 
2wf^2 + Mfg = i, 2w\H2 + 4u;f8 = ^i-4 = 0, 
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(16) 



There remains only three independent parameters in 
this expression, wf.g, W5-81 and wf_^. In the follow sec- 
tion we choose these unknowns in order to minimise the 
spurious velocity contribution. 



III. DETERMINING A UNIQUE EQUILIBRIUM 
(14) DISTRIBUTION 



In this section, we explicitly calculate the spurious 
force per unit volume G (see Eq. (IT6l) ) for the case of 
a liquid drop of radius R. If we take the origin to lie 
at the centre of the drop, then the density p is solely a 
function of distance from that origin r — \/x'^ + y^. Tay- 
lor expanding the po stencil (see Eq. ([TU])) we find the 
contribution to the force from this term is given by 

= - (dx + lAx^dl + 2wlsAx^dxdl) PO. (17) 

Transforming from Cartesian into polar coordinates is 
achieved using the relations 



dx 



-d d 

Ur , Uy 



(18) 



By sequentially substituting these operators and per- 
forming derivatives, Eq. (|17p can be rewritten 



Gl 



-xDrPo 



'1„3 
^6 



(i + 2<g) Ax^DIpo 
xy'') Ax'dIpq, 



2<8 



(19) 



where we define Dr — ^dr. By symmetry, the y compo- 
nent can be obtained by interchanging the x and y labels 
in this expression. The force can be decomposed into 
two terms; a term parallel and a term perpendicular to 
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the interface. The perpendicular contribution results in a 
small deviation in the Laplace pressure difference across 
the interface. The parallel term cannot be corrected for 
in this way and thus it is responsible for inducing spuri- 
ous flows. 

A parallel unit vector is given by n|| = (— y, x) and, 
therefore, this tangential contribution can be calculated 
using 



{x^y - xy^) (i - 2wl^) D^po. 



(20) 



This is zero provided that w^_g — . An analysis of other 
terms in Eq. (jl6p can be performed in a similar way. For 
instance, the force contribution from the Laplacian term 
is given by 



(9. 



2wl,Ax^d,d^y) 



^{dt + d^^)p + DAx'dld'^p) 



(21) 



By repeating the process that was used to derive Eq. 
P0|) . we find that this contribution vanishes provided 
that Wg g = j2 and D = ^. 

When transformed into polar coordinates, the tangen- 
tial force from the Maf3 terms in Eq. (fT6|) is given by 



{x^y 



xy 



+ [x^y - xy^) (- 



12 



12 



2<8 



Ax'^Dl{Drpf 
Ax^Dl{Drpf 



+ (i^y - xy'') (i - 2B) Ax^D,. [{Drp){D^^p)] 
+ {x^y - xy^) (i - 2B) Ax^{Drp)iDlp)\ . (22) 

The last two lines on the right hand side become zero 
when B = j^. Generally, it is not possible to make the 
first two lines simultaneously zero. However, it turns out 
that the first term dominates over the second, and so the 



best choice is 



The reason for this is that the 



^5-8 ~ 24- 

width of the interface is much smaller than the radius of 
curvature. The density p is approximately constant in 
the bulk regions but varies sharply in the interface. If 
we denote the width of the interface to be W then the 
largest value for a derivative can be typically obtained 
using ~ -(^F. Since the operator Dr appears one more 
time on the first line than the second, and it contains 
an extra factor of x^ or y^, then we expect the ratio in 
the magnitude of the first two lines to be approximately 

~ ^ W' -'^^ f&ct, a detailed analysis explicitly cal- 
culating the two functions based on a hyperbolic tangent 
interface profile reveals that their maxima differ by a fac- 
tor 3-^. Thus provided W the second line will be 
negligible compared to the first. 

Now that we have obtained a unique choice for the 
equilibrium, it is interesting to note that, to the best 
of our knowledge, none of the previously proposed free 
energy lattice Boltzmann schemes make this optimum 
choice. For example, Inamuro et al. [5] choose w^.g, — 
and Desplat et al. [8] choose w^.s = ~j2- 



;(b)i 



FIG. 1: The steady state velocity profile around a droplet 
using (a) a standard choice of equilibrium distribution and 
(b) an improved choice of equilibrium. 



IV. NUMERICAL RESULTS 

To test the predictions made in the previous section, we 
perform simulations on a grid of size 100 x 100. Param- 
eters used were a = 6 = and T = 0.56, leading 
to liquid-gas phase separation with densities pi = 4.54 
and pg — 2.57. The interfacial tension was set using 
K, = 0.025, giving an interface width of approximately 3 
lattice sites. 

A drop of radius R — 2b was initialised at the centre of 
the system and simulations were run for 10'* time-steps to 
allow steady state conditions to be reached. Figure [1]^ a) 
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1 
Parameter Value (x1/12) 



with wf.g. It clearly reaches a minimum very close to 
that predicted theoretically (wf.g = j^)- The spurious 
velocity never reaches exactly zero because our analysis 
only considered terms up to 0(9^) in the Taylor series 
expansion for the stencils (Eqs. (|10p and (fTTj) ). In real- 
ity, higher order terms also induce spurious velocities but 
these terms will be ^ smaller, and so have much less 
effect provided the interface width is reasonably large. 

The other curves in Fig. [2] show minima which cor- 
respond well with the values w\_^ — j^, wfjl — — 
B — j^, and D — ^ predicted in section Hill Note that 
when the simplest choice for calculating the derivatives is 
used (see Eq. ([9|)), the spurious velocities are ~ 10 times 
larger than for the optimum choice (this corresponds to 
B = in Fig. HJb)). 




V. WHAT HAPPENS WHEN r / 1 



1 2 3 

Parameter Value (x1/12) 



FIG. 2: The maximum spurious velocity as a function of (a) 
the equilibrium distribution parameters Wg,^ (solid line) , wl_g 
(dashed line), uif,^ (dotted line), and (b) the stencil param- 
eters B (solid line), D (dashed line), and the force stencil 
parameter F (dotted line). 



shows the flow profile around the drop for a typical set 
of parameters. We clearly observe eight vortices in the 
gas phase surrounding the curved interface of the drop. 
Fig. [ijb) shows the dramatic reduction in the spurious 
flow when the best parameter choice is used. 

To verify that the we have, indeed, obtained an op- 
timum choice of parameters, we show that the spurious 
velocity is minimised for each of the parameters sepa- 
rately. The effect of changing one parameter in isola- 
tion was found numerically by fixing all other degrees of 
freedom and scanning the chosen parameter's value over 
some range. This scanning procedure was performed suf- 
ficiently slowly to be quasi-static. Figure [2] shows the re- 
sults. The spurious velocity on the y-axis is defined to 
be the maximum velocity magnitude in the system. The 
solid curve in Fig. [DJa) shows how this velocity varies 



To obtain Eq. ([16)) we assumed that r = 1 and so 
the lattice Boltzmann equation reduced to /^(r, t + At) = 
/j^''(r — e^At, t). The more general case can be calculated 
under steady state conditions by sequentially substitut- 
ing Eq. (|4]) back into the fi term on the right hand side. 
This gives 



/r(r-e,Ai) + /r(r-2e,At)(l 



-H/r(r - 3e,Ai) (l 



(23) 



Therefore, /i(r) can be expressed in terms of the equi- 
librium distributions along lines of points radiating out 
following the velocity vector directions. The magnitude 
of these contributions decrease by a factor z — {1 — ^) 
for each step away. The stencils in Eq. fTB]) are no longer 
finite in size. For instance, the inner 5x5 region of the 
po stencil now looks like 
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(24) 



Converting this into continuous operators gives 



where the sum S is 



i=l 

- + 6r^ 



IPo,(25) 



(26) 



Since S is simply a numerical factor multiplying all the 
0{d'^) terms then it will also pre- multiply the spurious 
force expressions in Eqs. (P(7| and (P^ . Such a change 
does not alter the optimum choice of equilibrium when 
r^l. 
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FIG. 3: The variation in spurious velocity as a function of 
T for the standard LB (sohd line) and for the new forcing 
method (dashed line). 



The solid line in Fig. [3] shows how the numerically 
calculated spurious velocities depend on r using the op- 
timum choice for all other parameters. The function S 
passes through zero when t — ^ + = 0.789. This 
condition was calculated by Swift et al. [l[ using a dif- 
ferent method. It does not correspond exactly with the 
minimum of the curve because higher order spurious ve- 
locities become important in this region. 



As r is increased the spurious velocities rapidly in- 
crease in magnitude. These are principally generated by 
the small term on the second line of Eq. being mul- 
tiplied by the very large numerical factor S, which grows 



as T" 



VI. USING A FORCING METHOD 



Rather than incorporate the problematic Mafi terms 
into the equilibrium distribution, it is also possible to 
put them into a body force. The term Fi in Eq. Q is 
given by 



from the equilibrium distribution ^ and replaced by 



= -2^xiMxx - Myy) - dyM^y 
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(28) 



in the x direction, gy may be obtained by interchanging 
the labels x and y and transposing the stencils. Such a 
procedure leaves the continuum Navier-Stokes equation 
unchanged. 

This method has the advantage of allowing extra de- 
grees of freedom in choosing the stencils as compared to 
the standard lattice Boltzmann. In particular, we can 
choose to have a symmetry between the derivatives in 
the X and y directions {i.e. the y stencil can be ob- 
tained by transposing the x). By comparison, Eq. 
clearly cannot have this property. This improvement in 
the isotropy of the governing equation helps to reduce 
spurious velocities further. 

The dotted line in Fig. [Ifb) shows numerical results 
of how the spurious velocities change as a function of the 
stencil parameter F. The minimum of this curve lies at 
F = corresponding to a standard choice. By compar- 
ing the magnitude of the spurious velocity at this point 
with the minima from the other curves, we conclude that 
the forcing method leads to a further ~ 5 fold reduction, 
giving a typical value of ~ 2 x 10~^c. 

Another advantage of using forcing is shown in Fig. [S] 
As T is increased the spurious velocities normally become 
non-negligible due to the large numerical factor S in Eq. 
(|26p multiplying the otherwise small contribution from 
the second line in Eq. P^ . In the forcing method this 
term goes to zero allowing for accurate simulation of more 
viscous systems. 

In general, the disadvantages of using forcing meth- 
ods are that they make boundary conditions more com- 
plicated and, if being run on a parallel computer, re- 
quire more information to be passed between computer 
micro-processors. This is because the standard two di- 
mensional lattice Boltzmann method only requires infor- 
mation from the surrounding 8 points to update each 
lattice site, whereas the forcing method requires infor- 
mation from 24 points. 



F 



3 
2c^ 



(27) 



where g is a body force that now appears on the right 
hand side of the Navier-Stokes equation (|A7|1 . In this new 



forcing scheme, the wf^, wf , and w^ terms are removed 



VII. EXTENSION TO 3D LATTICE 
BOLTZMANN SCHEMES 

A number of different lattice Boltzmann schemes have 
been proposed for simulating 3D systems using 15, 19 
or 27 lattice velocities. In this paper we find that 19 
lattice vectors are necessary to ensure the reduction in 
spurious velocities. One way to define the velocity vectors 
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in this model is the following: ei_6 lie along the nearest 
neighbour directions 




c -c 
c -c 
c -c 



and 67-18 are in the 12 square diagonal directions 




c -c c -c c -c c -c 
c +c -c -c c -c c -c 
c c -c -c c c -c -c 



Analogous to the definitions for the gradient and 
Laplacian stencils given in Eqns. (fTU|) and (dJ), we define 
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(31) 



Note that for a system one lattice unit wide this equilib- 
rium reduces to the 2D result. 





£> 
BCD 
D 



(29) 



where E = -6C-12D, 2A+8B = 1, C + 4D = 1 and the 
left, middle, and right matrices show slices of the stencil 
when Czi — c,0, and — c, respectively. 

In three dimensions, additional terms containing wf^, 
wf^, and wf^ appear in the equilibrium distribution 
By using the same procedure as in section IIIIl this dis- 
tribution can be uniquely defined. One additional com- 
plication in the three dimensional case is that there is no 
longer a single vector defining a tangent to the surface 
of a drop. Instead we use three vectors ny — {—y,x,0), 
n| — (0, ~z,y), and n| — {—z,0,x) and require that the 
spurious force parallel to each one of these is zero. For 
instance, the previous expression in Eq. (|22|) becomes 
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14) 



+ (x^'y - xy'') (i - 2B) Ax" Dr{Drp){Dlp) 



+ {yx^z^ - xy^z^) (i - 2B) Ax" Dr{Drp){Dlp) 



(30) 



For the first line to be zero then 



W7 



= — Jj. The second 



^7-10 — 24- 

line, therefore, is zero only when 'wff_i^ = j^. The last 
two lines vanish when B = j^. 

A summary of all parameters obtained using this pro- 



VIII. SUMMARY AND CONCLUSIONS 

In this paper we analysed the spurious velocities from 
two different methods: a standard lattice Boltzmann 
scheme and a new forcing method. 

Firstly, we calculated the spurious forces which origi- 
nate when the continuous operators in the Navier-Stokes 
equation are replaced by stencils (in other word the con- 
tribution from the next order in the Taylor series expan- 
sion of the stencils). Secondly, we identify that spurious 
velocities result from the component of these spurious 
forces acting parallel to the interface. Finally, we find 
that by making a suitable choice of the equilibrium dis- 
tribution and stencils we were able to set these parallel 
forces to zero (up to fourth order in the derivatives). 

In 2D, the best choice of stencils for calculating the 
derivatives and the Laplacian are: 
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Using the standard lattice Boltzmann model the equi- 
librium is given by Eq. ([5]) , where the optimum choice of 
parameters is 7«i_4 



p 



<8 



1 11 1 ,xx 

12' "^5-8 
VV ^ _ 
1-2 



W5-8 



P 



5-8 — 

= = = wf4 = i 

3_4 - wi:2 = -h "^1-4 = 0' aiid ^5-8 = 3- In 3D 
the corresponding results are summarised in Eqns. (j29|) 
and (|5T|). 

One way to improve spurious velocities further is to 
remove the Map = ndapdpp terms from the equilibrium 
distribution and implement them as a body force. This 
force is then explicitly calculated by taking derivatives 
of Map using the stencil in Eq. ^ (or Eq. ^ in 
the 3D case). The additional symmetry in the resulting 
equations leads to a further reduction in spurious velocity 
size. 
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APPENDIX A: THE MOMENTS 

To conserve mass and momentum the first two con- 
straints on the equilibrium distribution must be 



(Al) 



The higher order moments of f^'^ are chosen such that the 
resulting continuum equations describe the dynamics of 
a non-ideal fluid. A suitable choice is 

ffeiaeip = Pap + pUaUp + A (uaBpp 

i 

+U0daP + Sapu^d^p) , (A2) 



i 

where 



will become the shear and bulk kinematic viscosities, re- 
spectively. The speed of sound is given by = 
where po is the fluid pressure The term involving A 



on the right hand side of Eq. (jA2p is necessary to ensure 
Galilean invariance For an ideal gas with = ^ it 
is zero, but in the more general case it must be included. 
The moments of the forcing term Fi are defined by 

i i 

FiEiaeifj = Ai (uQ.g/3 -I- upga) , (A5) 

i 

where g is the body force per unit volume acting on the 
fluid. 

By applying the Chapman-Enskog expansion to the 
lattice Boltzmann equation ([I]) we obtain the conti- 
nuity equation for the total density 

dtp + da{pva) = Q, (A6) 
and the Navier-Stokes equation for the fluid momentum 

+df3 {vp {dpVa + daVp) + XpSapd^Vj) , (A7) 



where the fluid velocity is defined by 



v = u+f g. 



(A8) 



Note that this definition differs slightly from the lattice 
fluid velocity ([3|) in the case when the body force is non- 
zero. It is V, and not u, which is used to calculate the 
spurious velocities in section IVTl 
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